<!DOCTYPE html>
<html>

<head>
  <meta charset="utf-8" />
  <meta name="viewport" content="width=device-width, initial-scale=1">
  <title>Exploring Bayesian Optimization</title>
  <script defer src="js/template.v2.js"></script>
  <link rel="stylesheet" type="text/css" href="css/styles.css">
  <script defer src="js/hider.js"></script>
  <script defer src="js/gif-slider.js"></script>

</head>

<body onload="appendInputButtons();">

  <d-front-matter>
    <script id='distill-front-matter' type="text/json">
      {
        "title": "Exploring Bayesian Optimization",
        "description": "How to tune hyperparameters for your machine learning model using Bayesian optimization.",
        "authors": [{
            "author": "Apoorv Agnihotri",
            "authorURL": "https://apoorvagnihotri.github.io/",
            "affiliations": [{
              "name": "Indian Insitute of Technology Gandhinagar",
              "affiliationURL": "https://www.iitgn.ac.in/"
            }]
          },
          {
            "author": "Nipun Batra",
            "authorURL": "https://nipunbatra.github.io/",
            "affiliations": [{
              "name": "Indian Insitute of Technology Gandhinagar",
              "affiliationURL": "https://www.iitgn.ac.in/"
            }]
          }
        ],
        "katex": {
          "delimiters": [{
            "left": "$$",
            "right": "$$",
            "display": false
          }]
        }
      }
    </script>
  </d-front-matter>

  <d-title style="padding-bottom: 0">
    <p>Breaking Bayesian Optimization into small, sizeable chunks.</p>
  </d-title>

  <d-byline></d-byline>

  <d-article style="overflow-x: unset;">
    <p>
      Many modern machine learning algorithms have a large number of hyperparameters. To effectively use these algorithms, we need to pick good hyperparameter values.

      In this article, we talk about Bayesian Optimization, a suite of techniques often used to tune hyperparameters. More generally, Bayesian Optimization can be used to optimize any black-box function.
    </p>
    <p>

    </p>

    <!-- <p>
      We first study using Bayesian Optimization to maximize (optimize) a black-box function.
    </p> -->

    <h1>Mining Gold!</h1>
    <p>
      Let us start with the example of gold mining. Our goal is to mine for gold in an unknown land<d-footnote>Interestingly, our example is similar to one of the first use of Gaussian Processes (also called kriging)<d-cite key="goldKridge"></d-cite>, where Prof. Krige modeled the gold concentrations using a Gaussian Process.</d-footnote>.
      For now, we assume that the gold is distributed about a line. We want to find the location along this line with the maximum gold while only drilling a few times (as drilling is expensive).
    </p>

    <p>
      Let us suppose that the gold distribution <d-math>f(x)</d-math> looks something like the function below. It is bi-modal, with a maximum value around <d-math>x = 5</d-math>. For now, let us not worry about the X-axis or the Y-axis units.
    </p>

    <figure class="smaller-img">
      <d-figure><img src="images/MAB_gifs/GT.svg" /></d-figure>
    </figure>

    <p>
      Initially, we have no idea about the gold distribution. We can learn the gold distribution by drilling at different locations. However, this drilling is costly. Thus, we want to <!-- <strong> -->minimize the number of drillings required<!-- </strong> --> while still <!-- <strong> -->finding the location of maximum gold quickly<!-- </strong> -->.
    </p>

    <p>
      We now discuss two common objectives for the gold mining problem.
    </p>

    <ul>
      <li>
        <p>
          <strong>Problem 1: Best Estimate of Gold Distribution (Active Learning)</strong><br/>
          In this problem, we want to accurately estimate the gold distribution on the new land.<!--  along the line, using a small number of drillings.  --> We can not drill at every location due to the prohibitive cost. Instead, we should drill at locations providing <strong>high information</strong> about the gold distribution. This problem is akin to
          <strong>
            Active Learning<d-cite key="settles2009active,Tong2001"></d-cite>
          </strong>.
        </p>
      </li>

      <li>
        <p>
          <strong>Problem 2: Location of Maximum Gold (Bayesian Optimization)</strong><br/>
          In this problem, we want to find the location of the maximum gold content. We, again, can not drill at every location. Instead, we should drill at locations showing <strong>high promise</strong> about the gold content.<!-- , using a small number of drillings. --> This problem is akin to
          <strong>
            Bayesian Optimization<d-cite key="humanOut,nandoBOtut"></d-cite>
          </strong>.
        </p>
      </li>
    </ul>

    <p>
      We will soon see how these two problems are related, but not the same.
    </p>



    <h2>Active Learning</h2>

    <p>
      For many machine learning problems, unlabeled data is readily available. However, labeling (or querying) is often expensive. As an example, for a speech-to-text task, the annotation requires expert(s) to label words and sentences manually. Similarly, in our gold mining problem, drilling (akin to labeling) is expensive. 
    </p>

    <p>
      Active learning minimizes labeling costs while maximizing modeling accuracy. While there are various methods in active learning literature, we look at <strong>uncertainty reduction</strong>. This method proposes labeling the point whose model uncertainty is the highest. Often, the variance acts as a measure of uncertainty.
    </p>

    <p>
      Since we only know the true value of our function at a few points, we need a <em>surrogate model</em> for the values our function takes elsewhere. This surrogate should be flexible enough to model the true function. Using a Gaussian Process (GP) is a common choice, both because of its flexibility and its ability to give us uncertainty estimates
      <d-footnote>
          Gaussian Process supports setting of priors by using specific kernels and mean functions. One might want to look at this excellent Distill article<d-cite key="görtler2019a"></d-cite> on Gaussian Processes<d-cite key="Rasmussen2004"></d-cite> to learn more.<br><br>
          Please find <a href="https://youtu.be/EnXxO3BAgYk">this</a> amazing video from Javier González on Gaussian Processes. 
      </d-footnote>.
    </p>

    <p>
      Our surrogate model starts with a prior of <d-math>f(x)</d-math> &#8212 in the case of gold, we pick a prior assuming that it’s smoothly distributed
      <d-footnote>
        Specifics: We use a Matern 5/2 kernel due to its property of favoring doubly differentiable functions. <!-- In contrast, Matern 3/2 favors singly differentiable functions. --> See <a href="http://www.gaussianprocess.org/gpml/">Rasmussen and Williams 2004</a> and <a href="https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.Matern.html">scikit-learn</a>, for details regarding the Matern kernel.
      </d-footnote>.
      As we evaluate points (drilling), we get more data for our surrogate to learn from, updating it according to Bayes’ rule.
    </p>

    <figure class="l-page">
      <d-figure><img src="images/MAB_gifs/prior2posterior.png" /></d-figure>
    <figcaption>
      Each new data point updates our surrogate model, moving it closer to the ground truth. The black line and the grey shaded region indicate the mean <d-math>(\mu)</d-math> and uncertainty <d-math>(\mu \pm \sigma)</d-math> in our gold distribution estimate before and after drilling.
    </figcaption>
    </figure>

    <p>
      In the above example, we started with uniform uncertainty. But after our first update, the posterior is certain near <d-math>x = 0.5</d-math> and uncertain away from it. We could just keep adding more training points and obtain a more certain estimate of <d-math>f(x)</d-math>.
    </p>

    <p>
      However, we want to minimize the number of evaluations. Thus, we should choose the next query point "smartly" using active learning. Although there are many ways to pick smart points, we will be picking the most uncertain one.
    </p>

    <p>
      This gives us the following procedure for Active Learning:
    </p>
  <!-- 
    <h3>Surrogate Model</h3>
    <p>
      Active learning (along with Bayesian Optimization, which we see later) employs a surrogate model for modeling the unknown true function <d-math>f(x)</d-math>. The surrogate ideally models the true function closely. 
    </p>

    <h3>Bayesian Update</h3>
    <p>
      Every evaluation (drilling) of <d-math>f(x)</d-math> gives the surrogate model more data to learn. We can apply Bayes rule to obtain the surrogate posterior. At the end of an evaluation, the posterior becomes the prior for the next evaluation.
    </p>

    <p>
      Gaussian Processes are commonly used as surrogate models. One can set priors of a Gaussian Process (GP) by using specific kernels and mean functions. Moreover, GPs provide uncertainty estimates, which we leverage in both active learning and Bayesian Optimization.
    </p>

    <h3>Gaussian Processes</h3>

    <p>
      One might want to look at this excellent Distill article<d-cite key="görtler2019a"></d-cite> on Gaussian Processes<d-cite key="Rasmussen2004"></d-cite>.
      We use Gaussian Process regression to model the gold distribution.
    </p>

    <p>
      Let us visualize our true function <d-math>f(x)</d-math>. The gold distribution is bi-modal, with a maximum value around <d-math>x = 5</d-math>. For now, let us not worry about the X-axis or the Y-axis units.
    </p>
    <figure class="smaller-img">
      <d-figure><img src="images/MAB_gifs/GT.svg" /></d-figure>
    </figure>


    <h4 id="priormodel">Prior Model</h4>

    <p>
      We define a prior over a set of functions based on our initial beliefs of the black-box function<d-cite key="Rasmussen2004"></d-cite>. We consider the gold distribution to be smooth by setting the  prior appropriately. <d-footnote>We use a Matern 5/2 kernel due to its property of favoring doubly differentiable functions. In contrast, Matern 3/2 favors singly differentiable functions.

      See <a href="http://www.gaussianprocess.org/gpml/">Rasmussen and Williams 2004</a> and <a href="https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.Matern.html">scikit-learn</a>, for details regarding the Matern kernel.
    </d-footnote> The black line and the grey shaded region indicate the mean (<d-math>\mu</d-math>) and uncertainty  (<d-math>\mu \pm \sigma </d-math>) in our gold distribution estimate before drilling.

    </p>

    <figure>
      <d-figure><img src="images/MAB_gifs/prior.svg" /></d-figure>
    </figure>

    <h4 id="addingtrainingdata">Adding Training Data</h4>

    <p>
      Let us now add the point <d-math>< x = 0.5, \ y = f(0.5) ></d-math> to the training set.
    </p>

    <figure>
      <d-figure><img src="images/MAB_gifs/posterior.svg" /></d-figure>
    </figure>

    <p>
      We see our surrogate's posterior is certainty near <d-math>x = 0.5</d-math> and uncertain away from it. We can similarly add more training points and obtain a more certain estimate of <d-math>f(x)</d-math>. However, we want to minimize the number of evaluations. Thus, we should choose the next query point "smartly" using active learning. !-- <d-footnote>We can transform the data and fit a GP over <d-math>\log\left(f(x)\right)</d-math> instead of <d-math>f(x)</d-math> to ensure the  predictions are non-negative</d-footnote> .  --
    We now discuss the key idea of active learning.</p>
    <h3 id="activelearningprocedure">Active Learning Procedure
    </h3>
  -->

    <ol>
      <li>Choose and add the point with the highest uncertainty to the training set (by querying/labeling that point)</li>
      <li>Train on the new training set</li>
      <li>Go to #1 till convergence or budget elapsed</li>
    </ol>

    <p>
      Let us now visualize this process and see how our posterior changes at every iteration (after each drilling).
    </p>
    <figure class="gif-slider">
      <d-figure><img src="images/active-gp-img/0.png"/></d-figure>
    </figure>


    <p>The visualization shows that one can estimate the true distribution in a few iterations. Furthermore, the most uncertain positions are often the farthest points from the current evaluation points. At every iteration, active learning <strong>explores</strong> the domain to make the estimates better.
    </p>

    <h2 id="bayesianoptimization">Bayesian Optimization</h2>
    <p>
      In the previous section, we picked points in order to determine an accurate model of the gold content. But what if our goal is simply to find the location of maximum gold content? Of course, we could do active learning to estimate the  true function accurately and then find its maximum. But that seems pretty wasteful &#8212 why should we use evaluations improving our estimates of regions where the function expects low gold content when we only care about the maximum?
    </p>

    <P>
      This is the core question in Bayesian Optimization: "Based on what we know so far, which point should we evaluate next?" Remember that evaluating each point is expensive, so we want to pick carefully! In the active learning case, we picked the most uncertain point, exploring the function. But in Bayesian Optimization, we need to balance exploring uncertain regions, which might unexpectedly have high gold content, against focusing on regions we already know have higher gold content (a kind of exploitation).
    </P>

    <p>
      We make this decision with something called an acquisition function. Acquisition functions are heuristics for how desirable it is to evaluate a point, based on our present model<d-footnote>More details on acquisition functions can be accessed at on this <a href="https://botorch.org/docs/acquisition">link</a>.</d-footnote>. We will spend much of this section going through different options for acquisition functions.
    </p>

    <p>
      This brings us to how Bayesian Optimization works. At every step, we determine what the best point to evaluate next is according to the acquisition function by optimizing it. We then update our model and repeat this process to determine the next point to evaluate.
    </p>

    <p>
      You may be wondering what’s “Bayesian” about Bayesian Optimization if we’re just optimizing these acquisition functions. Well, at every step we maintain a model describing our estimates and uncertainty at each point, which we update according to Bayes’ rule<d-cite key="nandoBOLoop"></d-cite> at each step. Our acquisition functions are based on this model, and nothing would be possible without them!
    </p>
    
    <!-- 
      <p>
        In this problem, we aim to find the location of maximum gold content. One way to find the maximum would be first to run active learning to estimate the true function accurately, and then find its maximum. However, should we waste evaluations to improve the estimates, when we are only concerned with finding the maximum? Assuming that our black-box function is smooth, it might be a good idea to evaluate at or near locations where our surrogate model's prediction is the highest, i.e., to <strong>exploit</strong>. However, due to the limited evaluations, our model's predictions are inaccurate. One can improve the model by evaluating at points with high variance or performing <strong>exploration</strong>. Bayesian Optimization combines <strong>exploitation</strong> and <strong>exploration</strong>, whereas active learning solely <strong>explores</strong>.
      </p> 
    -->
    
    <h3>Formalizing Bayesian Optimization</h3>

      <span style="padding-bottom: 1em;">
        Let us now formally introduce Bayesian Optimization. Our goal is to find the location (<d-math>{x \in \mathbb{R}^d}</d-math>) corresponding to the global maximum (or minimum) of a function <d-math>f: \mathbb{R}^d \mapsto \mathbb{R}</d-math>.
        We present the general constraints in Bayesian Optimization and contrast them with the constraints in our gold mining example<d-footnote>The section below is based on the slides/talk from Peter Fraizer at Uber on Bayesian Optimization:
          <ul class="footnote-ul">
            <li> <a href="https://www.youtube.com/watch?v=c4KKvyWW_Xk">Youtube talk</a>,</li>
            <li> <a href="https://people.orie.cornell.edu/pfrazier/Presentations/2018.11.INFORMS.tutorial.pdf">slide deck</a></li>
          </ul>
        </d-footnote>.</span>

      <table>
        <tr>
          <th><h3>General Constraints</h3></th>
          <th><h3>Constraints in Gold Mining example</h3></th>
        </tr>
        <tr>
          <td><d-math>f</d-math>’s feasible set <d-math>A</d-math> is simple,
          e.g., box constraints.</td>
          <td>Our domain in the gold mining problem is a single-dimensional box constraint: <d-math>0 \leq x \leq 6</d-math>.</td>
        </tr>
        <tr>
          <td><d-math>f</d-math> is continuous but lacks special structure,
          e.g., concavity, that would make it easy to optimize.</td>
          <td>Our true function is neither a convex nor a concave function, resulting in local optimums.</td>
        </tr>
        <tr>
          <td><d-math>f</d-math> is derivative-free:
          evaluations do not give gradient information.</td>
          <td>Our evaluation (by drilling) of the amount of gold content at a location did not give us any gradient information.</td>
        </tr>
        <tr>
          <td><d-math>f</d-math> is expensive to evaluate:
          the number of times we can evaluate it
          is severely limited.</td>
          <td>Drilling is costly.</td>
        </tr>
        <tr>
          <td><d-math>f</d-math> may be noisy. If noise is present, we will assume it is independent and normally distributed, with common but unknown variance.</td>
          <td>We assume noiseless measurements in our modeling (though, it is easy to incorporate normally distributed noise for GP regression).</td>
        </tr>
      </table>

      <p>
        To solve this problem, we will follow the following algorithm:
      </p>

      <p>
      <ol>
        <li>
          We first choose a surrogate model for modeling the true function <d-math>f</d-math> and define its <b>prior</b>.
        </li>
        <li>
          Given the set of <b>observations</b> (function evaluations), use Bayes rule to obtain the <b>posterior</b>.
        </li>
        <li>
          Use an acquisition function <d-math>\alpha(x)</d-math>, which is a function of the posterior, to decide the next sample point <d-math>x_t = \text{argmax}_x \alpha(x)</d-math>.
        </li>
        <li>
          Add newly sampled data to the set of <b>observations</b> and goto step #2 till convergence or budget elapses.
        </li>
      </ol>
    </p>

      <!-- <p>
        Our gold mining problem is suited to use Bayesian Optimization. Let us introduce some additional topics before we run to get the maximal gold!
      </p> -->

    <h3>Acquisition Functions</h3>

    <p>
      Acquisition functions are crucial to Bayesian Optimization, and there are a wide variety of options<d-footnote>
        Please find <a href="https://www.cse.wustl.edu/~garnett/cse515t/spring_2015/files/lecture_notes/12.pdf">these</a> slides from Washington University in St. Louis to know more about acquisition functions.
      </d-footnote>. In the following sections, we will go through a number of options, providing intuition and examples.
    </p>

    <!-- <p>
      The core question is Bayesian Optimization is "Based on what we know so far, which point should we evaluate next?" Remember that evaluating each point is expensive, so we want to pick carefully! In the active learning case, we picked the most uncertain point, but in Bayesian Optimization we don't waste precious evaluations on uncertain regions that are not promising.
    </p>

    <p>
      Acquisition functions are heuristics<d-footnote>https://botorch.org/docs/acquisition</d-footnote> for how desirable it is to evaluate a point, based on our present model. At every step, we determine what the best point to evaluate next is according to the acquisition function. We then update our model, and repeat this process to determine the next point to evaluate. 
      <d-footnote>
        Please find <a href="https://www.cse.wustl.edu/~garnett/cse515t/spring_2015/files/lecture_notes/12.pdf">these</a> slides from Washington University in St. Louis to know more about acquisition functions.
      </d-footnote>.
    </p> -->

    <!-- 
      <p>
        Acquisition functions are crucial to Bayesian Optimization, and there are a wide variety of options. In the following sections, we’ll go through a number of options, providing intuition and examples.
      </p>

      <p>
        Our original optimization problem, <d-math>x^* = \text{argmax}_{x \in A} f(x)</d-math> is hard because <d-math>f</d-math> is <b>expensive</b> to evaluate.
        The idea of Bayesian Optimization is to <b>transform</b> the original optimization into a <b>sequence</b> of easier <b>inexpensive</b> optimizations of functions called <b>acquisition functions</b> (<d-math>\alpha(x)</d-math>).
        Intuitively, acquisition functions are heuristics<d-footnote>https://botorch.org/docs/acquisition</d-footnote> that evaluate the utility of a point for maximizing <d-math>f(x)</d-math><d-footnote>Please find <a href="https://www.cse.wustl.edu/~garnett/cse515t/spring_2015/files/lecture_notes/12.pdf">these</a> slides from Washington University in St. Louis to know more and the following </d-footnote>.
        At each step, we <strong>optimize</strong> the <strong>acquisition function</strong> to determine the next point to sample/query.
      </p>
    -->

    <!-- <p>
      Let us re-wind and link the things discussed thus far, by noting the steps of Bayesian Optimization<d-footnote>Please find <a href="https://youtu.be/EnXxO3BAgYk">this</a> amazing video from Javier González at The Gaussian Process Summer School 2019.</d-footnote> and explicitly highlighting the "Bayesian" in Bayesian Optimization.
    </p>
    <p>
      <ol>
        <li>
          We first choose a surrogate model for modeling the true function <d-math>f</d-math> and define its <b>prior</b>.
        </li>
        <li>
          Given the set of <b>observations</b> (function evaluations), use Bayes rule to obtain the <b>posterior</b>.
        </li>
        <li>
          Use an acquisition function <d-math>\alpha(x)</d-math>, which is a function of the posterior, to decide the next sample point <d-math>x_t = \text{argmax}_x \alpha(x)</d-math>.
        </li>
        <li>
          Add newly sampled data to the set of <b>observations</b> and goto step #2 till convergence or budget elapses.
        </li>
      </ol>
    </p>

    <p>
      Thus, the "Bayesian" in Bayesian Optimization is sequentially refining our surrogate's posterior (and thus uncertainty) with each evaluation via Bayes rule<d-cite key="nandoBOLoop"></d-cite>.
    </p>
    <p>Let us now look at a few common acquisition functions.</p> -->

    <h4>Probability of Improvement (PI)</h4>

    <p>
      This acquisition function chooses the next query point as the one which has the highest <em>probability of improvement</em> over the current max <d-math>f(x^+)</d-math>. Mathematically, we write the selection of next point as follows, 
    </p>
    <div class="desktop">
      <d-math block>
        x_{t+1} = argmax(\alpha_{PI}(x)) = argmax(P(f(x) \geq (f(x^+) +\epsilon)))
      </d-math>
    </div>
    <div class="mobile_device_480px">
      <d-math block>
        \begin{aligned}
        x_{t+1} & = argmax(\alpha_{PI}(x))\\
        & = argmax(P(f(x) \geq (f(x^+) +\epsilon)))
        \end{aligned}
      </d-math>
    </div>
    <p>
      where, 
      <ul style="list-style-type: none;">
        <li><d-math>P(\cdot)</d-math> indicates probability</li>
        <li><d-math>\epsilon</d-math> is a small positive number</li>
        <li>And, <d-math> x^+ = \text{argmax}_{x_i \in x_{1:t}}f(x_i)</d-math> where <d-math>x_i</d-math> is the location queried at <d-math>i^{th}</d-math> time step.</li>
      </ul>
    </p>
    <p> 
      Looking closely, we are just finding the upper-tail probability (or the CDF) of the surrogate posterior. Moreover, if we are using a GP as a surrogate the expression above converts to,
    </p>
    <d-math block>x_{t+1} = argmax_x \Phi\left(\frac{\mu_t(x) - f(x^+) - \epsilon}{\sigma_t(x)}\right)</d-math>
    <p>
      where, 
      <ul style="list-style-type: none;">
        <li><d-math>\Phi(\cdot)</d-math> indicates the CDF</li>
      </ul>
    </p>

    <p>
      The visualization below shows the calculation of <d-math>\alpha_{PI}(x)</d-math>. The orange line represents the current max (plus an <d-math> \epsilon</d-math>) or <d-math> f(x^+) + \epsilon</d-math>. The violet region shows the probability density at each point. The grey regions show the probability density below the current max. The "area" of the violet region at each point represents the "probability of improvement over current maximum". The next point to evaluate via the PI criteria (shown in dashed blue line) is <d-math>x = 6</d-math>.
    </p>

    <figure>
      <d-figure><img src="images/MAB_gifs/density_pi.png" /></d-figure>
    </figure>

    <h5>Intuition behind <span style='text-transform: none'><d-math>\epsilon</d-math></span> in PI</h5>

    <p>
      PI uses <d-math>\epsilon</d-math> to strike a balance between exploration and exploitation. <!-- In the following plot, we visualize how changing the value of <d-math>\epsilon</d-math> affects the values of our acquisition function, resulting in different choices of points to evaluate. We selected two points as our candidate points <d-math>x \in \{1.0, 5.0\}</d-math>, and show their tail probabilities in shaded green and orange, respectively. -->
      Increasing <d-math>\epsilon</d-math> results in querying locations with a larger <d-math>\sigma</d-math> as their probability density is spread.
    </p>

    <p>
      Let us now see the PI acquisition function in action. We start with <d-math>\epsilon=0.075</d-math>.
    </p>

    <figure class="gif-slider">
      <d-figure><img src="images/MAB_pngs/PI/0.075/0.png" /></d-figure>
    </figure>
    <p>
      Looking at the graph above, we see that we reach the global maxima in a few iterations<d-footnote>Ties are broken randomly.</d-footnote>.
      Our surrogate possesses a large uncertainty in <d-math>x \in [2, 4]</d-math> in the first few iterations<d-footnote>The proportion of uncertainty is identified by the grey translucent area.</d-footnote>.
      The acquisition function initially <strong>exploits</strong> regions with a high promise<d-footnote>Points in the vicinity of current maxima</d-footnote>, which leads to high uncertainty in the region <d-math>x \in [2, 4]</d-math>. This observation also shows that we do not need to construct an accurate estimate of the black-box function to find its maximum.
    </p>

    <figure class="gif-slider">
      <d-figure><img src="images/MAB_pngs/PI/0.3/0.png" /></d-figure>
    </figure>
    <p>
      The visualization above shows that increasing <d-math>\epsilon</d-math> to 0.3, enables us to <strong>explore</strong> more. However, it seems that we are exploring more than required.
    </p>

    <p>
      What happens if we increase <d-math>\epsilon</d-math> a bit more?
    </p>

    <figure class="gif-slider">
      <d-figure><img src="images/MAB_pngs/PI/3/0.png" /></d-figure>
    </figure>
    <p>
      We see that we made things worse! Our model now uses <d-math>\epsilon = 3</d-math>, and we are unable to exploit when we land near the global maximum. Moreover, with high exploration, the setting becomes similar to active learning.
    </p>

    <p> Our quick experiments above help us conclude that <d-math>\epsilon</d-math> controls the degree of exploration in the PI acquisition function.

    <h4 id="expectedimprovementei">Expected Improvement (EI)</h4>

    <p>
      Probability of improvement only looked at <em>how likely</em> is an improvement, but, did not consider <em>how much</em> we can improve. The next criterion, called Expected Improvement (EI), does exactly that<d-footnote>A good introduction to the Expected Improvement acquisition function is by <a href="https://thuijskens.github.io/2016/12/29/bayesian-optimisation/" target="_blank">this post</a> by Thomas Huijskens and <a
          href="https://people.orie.cornell.edu/pfrazier/Presentations/2018.11.INFORMS.tutorial.pdf" target="_blank">these slides</a> by Peter Frazier</d-footnote>!
      The idea is fairly simple &#8212 choose the next query point as the one which has the highest expected improvement over the current max <d-math>f(x^+)</d-math>, where <d-math> x^+ = \text{argmax}_{x_i \in x_{1:t}}f(x_i)</d-math> and <d-math>x_i</d-math> is the location queried at <d-math>i^{th}</d-math> time step.
    </p>

    <p>
      In this acquisition function, <d-math>t + 1^{th}</d-math> query point, <d-math>x_{t+1}</d-math>, is selected according to the following equation.
    </p>
    <d-math block>
      x_{t+1} = argmin_x \mathbb{E} \left( ||h_{t+1}(x) - f(x^\star) || \ | \ \mathcal{D}_t \right)
    </d-math>
    <p>
      Where, <d-math>f</d-math> is the actual ground truth function, <d-math>h_{t+1}</d-math> is the posterior mean of the surrogate at <d-math>t+1^{th}</d-math> timestep, <d-math>\mathcal{D}_t</d-math> is the training data <d-math>\{(x_i,
        f(x_i))\} \ \forall x \in x_{1:t}</d-math> and <d-math>x^\star</d-math> is the actual position where <d-math>f</d-math> takes the maximum value.
    </p>

    <p>
      In essence, we are trying to select the point that minimizes the distance to the objective evaluated at the maximum. Unfortunately, we do not know the ground truth function, <d-math>f</d-math>. Mockus<d-cite key="mockusEI"></d-cite> proposed
      the following acquisition function to overcome the issue.
    </p>

    <div class="desktop">
      <d-math block>
      x_{t+1} = argmax_x \mathbb{E} \left( {max} \{ 0, \ h_{t+1}(x) - f(x^+) \} \ | \ \mathcal{D}_t \right)
      </d-math>
    </div>
    <div class="mobile_device_480px">
      <d-math block>
        \begin{aligned}
          x_{t+1} = \ & argmax_x \mathbb{E} \\
          & \left( {max} \{ 0, \ h_{t+1}(x) - f(x^+) \} \ | \ \mathcal{D}_t \right)
        \end{aligned}
      </d-math>
    </div>

    <p>
      where <d-math>f(x^+)</d-math> is the maximum value that has been encountered so far. This equation for GP surrogate is an analytical expression shown below.
    </p>

    <div class="desktop">
      <d-math block>
        EI(x)=
        \begin{cases}
        (\mu_t(x) - f(x^+) - \epsilon)\Phi(Z) + \sigma_t(x)\phi(Z), & \text{if}\ \sigma_t(x) > 0 \\
        0, & \text{if}\ \sigma_t(x) = 0
        \end{cases}
      </d-math>
    </div>
    <div class="mobile_device_480px">
      <d-math block>
        EI(x)= \begin{cases}
          [(\mu_t(x) - f(x^+) - \epsilon) & \ \sigma_t(x) > 0 \\
          \quad * \Phi(Z)] + \sigma_t(x)\phi(Z),\\
          0, & \ \sigma_t(x) = 0
        \end{cases}
      </d-math>
    </div>
    <d-math block>Z= \frac{\mu_t(x) - f(x^+) - \epsilon}{\sigma_t(x)}</d-math>
    <p>
      where <d-math>\Phi(\cdot)</d-math> indicates CDF and <d-math>\phi(\cdot)</d-math> indicates pdf.
    </p>
    <p>From the above expression, we can see that <em>Expected Improvement</em> will be high when: i) the expected value of <d-math>\mu_t(x) - f(x^+)</d-math> is high, or, ii) when the uncertainty <d-math>\sigma_t(x)</d-math> around a point is high.


    <p> Like the PI acquisition function, we can moderate the amount of exploration of the EI acquisition function by modifying <d-math>\epsilon</d-math>.
    </p>

    <figure class="gif-slider">
      <d-figure><img src="images/MAB_pngs/EI/0.01/0.png" /></d-figure>
    </figure>

    <p>
      For <d-math>\epsilon = 0.01</d-math> we come close to the global maxima in a few iterations.  </p>

    <p>
      We now increase <d-math>\epsilon</d-math> to explore more.
    </p>

    <figure class="gif-slider">
      <d-figure><img src="images/MAB_pngs/EI/0.3/0.png" /></d-figure>
    </figure>
    <p>
      As we expected, increasing the value to <d-math>\epsilon = 0.3</d-math> makes the acquisition function explore more. Compared to the earlier evaluations, we see less exploitation. We see that it evaluates only two points near the global maxima.
    </p>

    <p>
      Let us increase <d-math>\epsilon</d-math> even more.
    </p>

    <figure class="gif-slider">
      <d-figure><img src="images/MAB_pngs/EI/3/0.png" /></d-figure>
    </figure>
    <p>
      Is this better than before? It turns out a yes and a no; we explored too much at <d-math>\epsilon = 3</d-math> and quickly reached near the global maxima. But unfortunately, we did not exploit to get more gains near the global maxima.
    </p>

    <div class="collapsible"><h4>PI vs Ei</h4><div class="collapsible-indicator"></div></div>
    <div class="content">
      <p>
        We have seen two closely related methods, The <em>Probability of Improvement</em> and the <em>Expected Improvement</em>.
      </p>

      <figure class="smaller-img">
        <d-figure><img src="images/MAB_gifs/Ei_Pi_graph/0.svg" /></d-figure>
      </figure>
      <p>
        The scatter plot above shows the policies' acquisition functions evaluated on different points<d-footnote>Each dot is a point in the search space. Additionally, the training set used while making the plot only consists of a single observation <d-math>(0.5, f(0.5))</d-math></d-footnote>.
        We see that <d-math>\alpha_{EI}</d-math> and <d-math>\alpha_{PI}</d-math> reach a maximum of 0.3 and around 0.47, respectively. Choosing a point with low <d-math>\alpha_{PI}</d-math> and high <d-math>\alpha_{EI}</d-math> translates to high risk<d-footnote>Since "Probability of Improvement" is low</d-footnote> and high reward<d-footnote>Since "Expected Improvement" is high</d-footnote>.
        In case of multiple points having the same <d-math>\alpha_{EI}</d-math>, we should prioritize the point with lesser risk (higher <d-math>\alpha_{PI}</d-math>). Similarly, when the risk is same (same <d-math>\alpha_{PI}</d-math>), we should choose the point with greater reward (higher <d-math>\alpha_{EI}</d-math>).
      </p>
    </div>

    <h3 id="thompsonsampling">Thompson Sampling</h3>

    <p>
      Another common acquisition function is Thompson Sampling <d-cite key="thompson"></d-cite>. At every step, we sample a function from the surrogate's posterior and optimize it. For example, in the case of gold mining, we would sample a plausible distribution of the gold given the evidence and evaluate (drill) wherever it peaks.
    </p>

    <p>
      Below we have an image showing three sampled functions from the learned surrogate posterior for our gold mining problem. The training data constituted the point <d-math>x = 0.5</d-math> and the corresponding functional value.
    </p>

    <figure>
      <d-figure>
        <img src="images/MAB_gifs/thompson.svg" />
      </d-figure>
    </figure>

    <p>
      We can understand the intuition behind Thompson sampling by two observations:
      <ul>
        <li>
          <p>
            Locations with high uncertainty (<d-math> \sigma(x) </d-math>) will show a large variance in the functional values sampled from the surrogate posterior. Thus, there is a non-trivial probability that a sample can take high value in a highly uncertain region. Optimizing such samples can aid <strong>exploration</strong>.
          </p>
          <p>
            As an example, the three samples (sample #1, #2, #3) show a high variance close to <d-math>x=6</d-math>. Optimizing sample 3 will aid in exploration by evaluating <d-math>x=6</d-math>.
          </p>

        </li>
        <li>
          <p>
            The sampled functions must pass through the current max value, as there is no uncertainty at the evaluated locations. Thus, optimizing samples from the surrogate posterior will ensure <strong>exploiting</strong> behavior.
          </p>

          <p>
            As an example of this behavior, we see that all the sampled functions above pass through the current max at <d-math>x = 0.5</d-math>. If <d-math>x = 0.5</d-math> were close to the global maxima, then we would be able to <strong>exploit</strong> and choose a better maximum.
          </p>

        </li>
      </ul>



      <figure class="gif-slider">
        <d-figure><img src="images/MAB_pngs/Thompson/0.png" /></d-figure>
      </figure>

      <p>
The visualization above uses Thompson sampling for optimization. Again, we can reach the global optimum in relatively few iterations.
      </p>

      <h3>Random</h3>

      <p>
        We have been using intelligent acquisition functions until now.
        We can create a random acquisition function by sampling <d-math>x</d-math>
        randomly.      </p>

      <figure class="gif-slider">
        <d-figure><img src="images/MAB_pngs/Rand/0.png" /></d-figure>
      </figure>
	<p> The visualization above shows that the performance of the random acquisition function is not that bad! However, if our optimization was more complex (more dimensions), then the random acquisition might perform poorly.
</p>
      <h3>Summary of Acquisition Functions</h3> <p>
        Let us now summarize the core ideas associated with acquisition functions: i) they are heuristics for evaluating the utility of a point; ii) they are a function of the surrogate posterior; iii) they combine exploration and exploitation; and iv) they are inexpensive to evaluate.</p>

      <div class="collapsible"><h4>Other Acquisition Functions</h4> <div class="collapsible-indicator"></div></div>
      <div class="content">

        <p>We have seen various acquisition functions until now. One trivial way to come up with acquisition functions is to have a explore/exploit combination.
        <p>

        <h3> Upper Confidence Bound (UCB) </h3>
        <p>
          One such trivial acquisition function that combines the exploration/exploitation tradeoff is a linear combination of the mean and uncertainty of our surrogate model. The model mean signifies exploitation (of our model's knowledge) and model uncertainty signifies exploration (due to our model's lack of observations).
          <d-math block>\alpha(x) = \mu(x) + \lambda \times \sigma(x)</d-math>
        </p>

        <p>
          The intuition behind the UCB acquisition function is weighing of the  importance between the surrogate's mean  vs. the surrogate's uncertainty. The <d-math>\lambda</d-math> above is the hyperparameter that can control the preference between exploitation or exploration.
        </p>

        <p>
          We can further form acquisition functions by combining the existing acquisition functions though the physical interpretability of such combinations might not be so straightforward. One reason we might want to combine two methods is to overcome the limitations of the individual methods.
        </p>

        <h3>Probability of Improvement + <d-math>\lambda \ \times</d-math> Expected Improvement (EI-PI)</h3>

        <p>
          One such combination can be a linear combination of PI and EI.

          We know PI focuses on the probability of improvement, whereas EI focuses on the expected improvement. Such a combination could help in having a tradeoff between the two based on the value of <d-math>\lambda</d-math>.
        </p>

        <h3>Gaussian Process Upper Confidence Bound (GP-UCB)</h3>

        <p>
          Before talking about GP-UCB, let us quickly talk about <strong>regret</strong>. Imagine if the maximum gold was <d-math>a</d-math> units, and our optimization instead samples a location containing <d-math>b < a</d-math> units, then our regret is
              <d-math>a -
              b</d-math>. If we accumulate the regret over <d-math>n</d-math> iterations, we get what is called <strong> cumulative regret. </strong> <br />

          GP-UCB's<d-cite key="gpucb"></d-cite> formulation is given by:
        </p>

        <d-math block>
          \alpha_{GP-UCB}(x) = \mu_t(x) + \sqrt{\beta_t}\sigma_t(x)
        </d-math>
        <p>
          Where <d-math>t</d-math> is the timestep.
        </p>

        <p>
          Srinivas et. al.<d-cite key="gpucbBounds"></d-cite> developed a schedule for <d-math>\beta</d-math> that they theoretically demonstrate to minimize cumulative regret.
        </p>

      </div>
      
      <h3>Comparison</h3>

      <p>
        We now compare the performance of different acquisition functions on the gold mining problem<d-footnote>To know more about the difference between acquisition functions look at <a href="https://www.cs.ubc.ca/~nando/540-2013/lectures/l7.pdf">these</a> amazing
          slides from Nando De Freitas</d-footnote>. We have used the optimum hyperparameters for each acquisition function.

          We ran the random acquisition function several times with different seeds and plotted the mean gold sensed at every iteration.
      </p>

      <figure>
        <d-figure><img src="images/MAB_gifs/comp.svg" /></d-figure>
      </figure>

      <p>
        The <em>random</em> strategy is initially comparable to or better than other acquisition functions<d-footnote>UCB and GP-UCB have been mentioned in the collapsible</d-footnote>. However, the maximum gold sensed by <em>random</em> strategy grows slowly. In comparison, the other acquisition functions can find a good solution in a small number of iterations. In fact, most acquisition functions reach fairly close to the global maxima in as few as three iterations.
      </p>

      <h2>Hyperparameter Tuning</h2>

      <p>Before we talk about Bayesian optimization for hyperparameter tuning<d-cite key="Snoek2012,NIPS2011_4443,Bergstra"></d-cite>, we will quickly differentiate between hyperparameters and parameters: hyperparameters are set before learning and the parameters are learned from the data. To illustrate the difference, we take the example of Ridge regression.
      </p>
      <div class="desktop">
      <d-math block>
        \hat{\theta}_{ridge} = argmin_{\theta\ \in \ \mathbb{R}^p} \sum\limits_{i=1}^{n} \left(y_i - x_i^T\theta \right)^2 + \lambda \sum\limits_{j=1}^{p} \theta^2_j
      </d-math>
    </div>
    <div class="mobile_device_480px">
      <d-math block>
        \begin{aligned}
        \hat{\theta}_{ridge} = & argmin_{\theta\ \in \ \mathbb{R}^p} \sum\limits_{i=1}^{n} \left(y_i - x_i^T\theta \right)^2 \\
        & + \lambda \sum\limits_{j=1}^{p} \theta^2_j
        \end{aligned}
      </d-math>
    </div>
      <p>
        In Ridge regression, the weight matrix <d-math>\theta</d-math> is the parameter, and the regularization coefficient <d-math>\lambda \geq 0</d-math> is the hyperparameter. <br />
        If we solve the above regression problem via gradient descent optimization, we further introduce another optimization parameter, the learning rate <d-math>\alpha</d-math>.
      </p>

      <p>The most common use case of Bayesian Optimization is <em>hyperparameter tuning</em>: finding the best performing hyperparameters on machine learning models.</p>

      <p>When training a model is not expensive and time-consuming, we can do a grid search to find the optimum hyperparameters. However, grid search is not feasible if function evaluations are costly, as in the case of a large neural network that takes days to train. Further, grid search scales poorly in terms of the number of hyperparameters.

      <p>We turn to Bayesian Optimization to counter the expensive nature of evaluating our black-box function (accuracy).</p>

      <h3 id="example1">Example 1 -- Support Vector Machine (SVM)</h3>

      <p>In this example, we use an SVM to classify on sklearn's moons dataset and use Bayesian Optimization to optimize SVM hyperparameters.

      <p>
        <ul>
          <li>
            <d-math>\gamma</d-math> -- modifies the behavior of the SVM's kernel. Intuitively it is a measure of the influence of a single training example<d-footnote>StackOverflow <a
                href="https://stackoverflow.com/questions/35848210/support-vector-machine-what-are-c-gamma">answer</a> for intuition behind the hyperparameters.</d-footnote>.
          </li>
          <li>
            <d-math>C</d-math> -- modifies the slackness of the classification, the higher the <d-math>C</d-math> is, the more sensitive is SVM towards the noise.
          </li>
        </ul>
      </p>

<p> Let us have a look at the dataset now, which has two classes and two features.</p>

      <figure class="smaller-img">
        <d-figure><img src="images/MAB_gifs/moons.svg" /></d-figure>
      </figure>

      <p>
        Let us apply Bayesian Optimization to learn the best hyperparameters for this classification task<d-footnote> <strong>Note</strong>: the surface plots you see for the Ground Truth Accuracies below were calculated for each possible hyperparameter for showcasing purposes only. We do not have these values in real applications.
        </d-footnote>. The optimum values for &#60<d-math>C, \ \gamma</d-math>&#62 have been found via running grid search at high granularity.
      </p>

      <figure class="gif-slider">
        <d-figure><img src="images/MAB_pngs/PI3d/0.05/0.png" /></d-figure>
      </figure>

      <p>Above we see a slider showing the work of the <em>Probability of Improvement</em> acquisition function in finding the best hyperparameters.</p>

      <figure class="gif-slider">
        <d-figure><img src="images/MAB_pngs/EI3d/0.01/0.png" /></d-figure>
      </figure>

      <p>Above we see a slider showing the work of the <em>Expected Improvement</em> acquisition function in finding the best hyperparameters.</p>

      <h3 id="comparison">Comparison</h3>

      <p>
        Below is a plot that compares the different acquisition functions. We ran the <em>random</em> acquisition function several times to average out its results.
      </p>

      <figure>
        <d-figure><img src="images/MAB_gifs/comp3d.svg" /></d-figure>
      </figure>

    <p>
        All our acquisition beat the <em>random</em> acquisition function after seven iterations. We see the <em>random</em> method seemed to perform much better initially, but it could not reach the global optimum, whereas Bayesian Optimization was able to get fairly close. The initial subpar performance of Bayesian Optimization can be attributed to the initial exploration.
    </p>

 

      <div class="collapsible"><h4>Other Examples</h4><div class="collapsible-indicator"></div></div>
      <div class="content">
        <h3>Example 2 -- Random Forest</h3>

        <p>
          Using Bayesian Optimization in a Random Forest Classifier.<d-cite key="scikit"></d-cite>
        </p>

        <p>
          We will continue now to train a Random Forest on the moons dataset we had used previously to learn the Support Vector Machine model. The primary hyperparameters of Random Forests we would like to optimize our accuracy are the <strong> number</strong> of
          Decision Trees we would like to have, the <strong>maximum depth</strong> for each of those decision trees.
        </p>
        <p>
          The parameters of the Random Forest are the individual trained Decision Trees models.
        </p>
        <p>
          We will be again using Gaussian Processes with Matern kernel to estimate and predict the accuracy function over the two hyperparameters.
        </p>

        <figure class="gif-slider">
          <d-figure><img src="images/MAB_pngs/RFPI3d/0.05/0.png"></d-figure>
        </figure>

        <p>
          Above is a typical Bayesian Optimization run with the <em>Probability of Improvement</em> acquisition function.
        </p>

        <figure class="gif-slider">
          <d-figure><img src="images/MAB_pngs/RFEI3d/0.5/0.png"></d-figure>
        </figure>

        <p>Above we see a run showing the work of the <em>Expected Improvement</em> acquisition function in optimizing the hyperparameters.</p>

        <figure class="gif-slider">
          <d-figure><img src="images/MAB_pngs/RFGP_UCB3d/1-2/0.png"></d-figure>
        </figure>

        <p>
          Now using the <em>Gaussian Processes Upper Confidence Bound</em> acquisition function in optimizing the hyperparameters.
        </p>

        <figure class="gif-slider">
          <d-figure><img src="images/MAB_pngs/RFRand3d/1-2/0.png"></d-figure>
        </figure>

        <p>Let us now use the Random acquisition function.</p>

        <figure>
          <d-figure><img src="images/MAB_gifs/RFcomp3d.svg"></d-figure>
        </figure>

        <p>
          The optimization strategies seemed to struggle in this example. This can be attributed to the non-smooth ground truth. This shows that the effectiveness of Bayesian Optimization depends on the surrogate's efficiency to model the actual black-box function. It is interesting to notice that the Bayesian Optimization framework still beats the <em>random</em> strategy using various acquisition functions.
        </p>

        <h3>Example 3 -- Neural Networks</h3>
        <p>
          Let us take this example to get an idea of how to apply Bayesian Optimization to train neural networks. Here we will be using <code>scikit-optim</code>, which also provides us support for optimizing function with a search space of categorical, integral, and real variables. We will not be plotting the ground truth here, as it is extremely costly to do so. Below are some code snippets that show the ease of using Bayesian Optimization packages for hyperparameter tuning.
        </p>

        <p>
          The code initially declares a search space for the optimization problem. We limit the search space to be the following:
          <ul>
            <li>
              batch_size -- This hyperparameter sets the number of training examples to combine to find the gradients for a single step in gradient descent. </br> 
              Our search space for the possible batch sizes consists of integer values s.t. batch_size = <d-math>2^i \ \forall \ 2 \leq i \leq 7 \ \& \ i \in \mathbb{Z}</d-math>.
            </li>
            <li>
              learning rate -- This hyperparameter sets the stepsize with which we will perform gradient descent in the neural network. </br>
              We will be searching over all the real numbers in the range <d-math>[10^{-6}, \ 1]</d-math>.
            </li>
            <li>
              activation -- We will have one categorical variable, i.e. the activation to apply to our neural network layers. This variable can take on values in the set <d-math>\{ relu, \ sigmoid \}</d-math>.
            </li>
          </ul>
        </p>

        <d-code block language="python">
          log_batch_size = Integer(
            low=2,
            high=7,
            name='log_batch_size'
          )
            lr = Real(
            low=1e-6,
            high=1e0,
            prior='log-uniform',
            name='lr'
          )
          activation = Categorical(
            categories=['relu', 'sigmoid'],
            name='activation'
          )

          dimensions = [
            dim_num_batch_size_to_base,
            dim_learning_rate,
            dim_activation
          ]
        </d-code>

        <p>
          Now import <code>gp-minimize</code><d-footnote><strong>Note</strong>: One will need to negate the accuracy values as we are using the minimizer function from <code>scikit-optim</code>.</d-footnote> from <code>scikit-optim</code> to perform the optimization. Below we show calling the optimizer using <em>Expected Improvement</em>, but of course we can select from a number of other acquisition functions.
        </p>

        <d-code block language="python">
          # initial parameters (1st point)
          default_parameters = 
            [4, 1e-1, 'relu']

          # bayesian optimization
            search_result = gp_minimize(
            func=train,
            dimensions=dimensions,
            acq_func='EI', # Expctd Imprv.
            n_calls=11,
            x0=default_parameters
          )
        </d-code>

        <figure class="smaller-img">
          <d-figure><img src="images/MAB_gifs/conv.svg"></d-figure>
        </figure>

        <p>
          In the graph above the y-axis denotes the best accuracy till then, <d-math>\left( f(x^+) \right)</d-math> and the x-axis denotes the evaluation number.
        </p>

        <p>
          Looking at the above example, we can see that incorporating Bayesian Optimization is not difficult and can save a lot of time. Optimizing to get an accuracy of nearly one in around seven iterations is impressive!<d-footnote>The example above has been inspired by <a href="https://github.com/Hvass-Labs/TensorFlow-Tutorials/blob/master/19_Hyper-Parameters.ipynb">Hvass Laboratories' Tutorial Notebook</a> showcasing hyperparameter optimization in TensorFlow using <code>scikit-optim</code>.</d-footnote>
        </p>

        <p>
          Let us get the numbers into perspective. If we had run this optimization using a grid search, it would have taken around <d-math>(5 \times 2 \times 7)</d-math> iterations. Whereas Bayesian Optimization only took seven iterations. Each iteration took around fifteen minutes; this sets the time required for the grid search to complete around seventeen hours!
        </p>
      </div>

      <h1 id="conclusions">Conclusion and Summary</h1>

      <p>
        In this article, we looked at Bayesian Optimization for optimizing a black-box function. Bayesian Optimization is well suited when the function evaluations are expensive, making grid or exhaustive search impractical. We looked at the key components of Bayesian Optimization. First, we looked at the notion of using a surrogate function (with a prior over the space of objective functions) to model our black-box function. Next, we looked at the "Bayes" in Bayesian Optimization -- the function evaluations are used as data to obtain the surrogate posterior. We look at acquisition functions, which are functions of the surrogate posterior and are optimized sequentially. This new sequential optimization is in-expensive and thus of utility of us. We also looked at a few acquisition functions and showed how these different functions balance exploration and exploitation. Finally, we looked at some practical examples of Bayesian Optimization for optimizing hyper-parameters for machine learning models.
      </p>

      <p>
        We hope you had a good time reading the article and hope you are ready to <strong>exploit</strong> the power of Bayesian Optimization. In case you wish to <strong>explore</strong> more, please read the <a href="#FurtherReading">Further Reading</a> section below. We also provide our <a href="https://github.com/distillpub/post--bayesian-optimization">repository</a> to reproduce the entire article.  
      </p>

  </d-article>

  <d-appendix>
    <h3 id="embracebayesianoptimization">Embrace Bayesian Optimization</h2>

      <p>
        Having read all the way through, you might have been sold on the idea about the time you can save by asking Bayesian Optimizer to find the best hyperparameters for your fantastic model. There are a plethora of Bayesian Optimization libraries available. We have linked a few below. Do check them out.
      </p>

      <ul>
        <li><a href="https://scikit-optimize.github.io/">scikit-optimize</a>
          <d-footnote>Really nice tutorial showcasing hyperparameter optimization on a neural network available at this <a href="https://github.com/Hvass-Labs/TensorFlow-Tutorials/blob/master/19_Hyper-Parameters.ipynb">link</a>.
          </d-footnote>
        </li>

        <li><a href="https://app.sigopt.com/docs/overview/python">sigopt</a></li>

        <li><a href="http://hyperopt.github.io/hyperopt/">hyperopt</a></li>

        <li><a href="https://github.com/HIPS/Spearmint">spearmint</a></li>

        <li><a href="https://github.com/Yelp/MOE">MOE</a></li>

        <li><a href="https://botorch.org/tutorials/">BOTorch</a></li>
        <li><a href="https://github.com/SheffieldML/GPyOpt">GPyOpt</a></li>
        <li><a href="https://github.com/dragonfly/dragonfly">DragonFly</a></li>
      </ul>

      <h3 id="ack">Acknowledgments</h2>

      <p>
        This article was made possible with inputs from numerous people. Firstly, we would like to thank all the Distill reviewers for their punctilious and actionable feedback. These fantastic reviews immensely helped strengthen our article. We further express our gratitude towards the Distill Editors, who were extremely kind and helped us navigate various steps to publish our work. We would also like to thank <a href="https://sgarg87.github.io/">Dr. Sahil Garg</a> for his feedback on the flow of the article. We would like to acknowledge the help we received from <a href="http://initiatives.iitgn.ac.in/writingstudio/wp/">Writing Studio</a> to improve the script of our article. Lastly, we sincerely thank <a href="https://colah.github.io/">Christopher Olah</a>. His inputs, suggestions, multiple rounds of iterations made this article substantially better.  
      </p>

    <h3 id="FurtherReading">Further Reading</h3>
    <ol>
      <li>
        <p>Using gradient information when it is available.</p>
        <ul>
          <li>
            Suppose we have gradient information available, we should possibly try to use the information. This could result in a much faster approach to the global maxima. Please have a look at the paper by Wu, et al.<d-cite key="BOwtGD"></d-cite> to know more.
          </li>
        </ul>
      </li>
      <li>
        <p>
          To have a quick view of differences between Bayesian Optimization and Gradient Descent, one can look at <a href="https://stats.stackexchange.com/q/161936">this</a> amazing answer at StackOverflow.
        </p>
      </li>
      <li>
        <p>
          We talked about optimizing a black-box function here. If we are to perform over multiple objectives, how do these acquisition functions scale? There has been fantastic work in this domain too! We try to deal with these cases by having multi-objective acquisition functions. Have a look at <a href="https://gpflowopt.readthedocs.io/en/latest/notebooks/multiobjective.html">this excellent</a> notebook for an example using <code>gpflowopt</code>.
        </p>
      </li>
      <li>
        <p>
          One of the more interesting uses of hyperparameters optimization can be attributed to searching the space of neural network architecture for finding the architectures that give us maximal predictive performance. One might also want to consider nonobjective optimizations as some of the other objectives like memory consumption, model size, or inference time also matter in practical scenarios.
        </p>
      </li>
      <li>
        <p>
          When the datasets are extremely large, human experts tend to test hyperparameters on smaller subsets of the dataset and iteratively improve the accuracy for their models. There has been work in Bayesian Optimization, taking into account these approaches<d-cite key="hyperband,largeBO"></d-cite> when datasets are of such sizes.
        </p>
      </li>
      <li>
        <p>
          There also has been work on Bayesian Optimization, where one explores with a certain level of "safety", meaning the evaluated values should lie above a certain security threshold functional value<d-cite key="SafeExplore"></d-cite>. One toy example is the possible configurations for a flying robot to maximize its stability. If we tried a point with terrible stability, we might crash the robot, and therefore we would like to explore the configuration space more diligently.
        </p>
      </li>
      <li>
        <p>
          We have been using GP in our Bayesian Optimization for getting predictions, but we can have any other predictor or mean and variance in our Bayesian Optimization.
        </p>
        <ul>
          <li>
            <p>
              One can look at <a href="http://aad.informatik.uni-freiburg.de/~hutter/ML3.pdf">this</a> slide deck by Frank Hutter discussing some limitations of a GP-based Bayesian Optimization over a Random Forest based Bayesian Optimization.
            </p>
          </li>
          <li>
            <p>
              There has been work on even using deep neural networks in Bayesian Optimization<d-cite key="NNbasedBO"></d-cite> for a more scalable approach compared to GP. The paper talks about how GP-based Bayesian Optimization scales cubically with the number of observations, compared to their novel method that scales linearly.
            </p>
          </li>
        </ul>
      </li>
      <li>
        <p>
          Things to take care when using Bayesian Optimization.
        </p>
        <ul>
          <li>
            <p>
              While working on the blog, we once scaled the accuracy from the range <d-math>[0, \ 1]</d-math> to <d-math>[0, \ 100]</d-math>. This change broke havoc as the Gaussian Processes we were using had certain hyperparameters, which needed
              to be scaled with the accuracy to maintain scale invariance. We wanted to point this out as it might be helpful for the readers who would like to start using on Bayesian Optimization.
            </p>
          </li>
          <li>
            <p>
              We need to take care while using Bayesian Optimization. Bayesian Optimization based on Gaussian Processes Regression is highly sensitive to the kernel used. For example, if you are using <a
                href="https://scikit-learn.org/stable/modules/generated/sklearn.gaussian_process.kernels.Matern.html">Matern</a> kernel, we are implicitly assuming that the function we are trying to optimize is first order differentiable.
            </p>
          </li>
          <li>
            <p>
              Searching for the hyperparameters, and the choice of the acquisition function to use in Bayesian Optimization are interesting problems in themselves. There has been amazing work done, looking at this problem. As mentioned previously in the post, there has
              been work done in strategies using multiple acquisition function<d-cite key="multiACQ"></d-cite> to deal with these interesting issues.
            </p>
          </li>
          <li>
            <p>
              A nice list of tips and tricks one should have a look at if you aim to use Bayesian Optimization in your workflow is from this fantastic post by Thomas on <a href="https://thuijskens.github.io/2016/12/29/bayesian-optimisation/">Bayesian
                Optimization with sklearn</a>.
            </p>
          </li>
        </ul>
      </li>
      <li>
        <p>
          Bayesian Optimization applications.
        </p>
        <ul>
          <li>
            <p>
              Bayesian Optimization has been applied to Optimal Sensor Set selection for predictive accuracy<d-cite key="sensorBO"></d-cite>.
            </p>
          </li>
          <li>
            <p>
              Peter Frazier in his <a href="https://www.youtube.com/watch?v=c4KKvyWW_Xk">talk</a> mentioned that Uber uses Bayesian Optimization for tuning algorithms via backtesting.
            </p>
          </li>
          <li>
            <p>Facebook<d-cite key="letham2019"></d-cite> uses Bayesian Optimization for A/B testing.
          </li>
          <li>
            <p>
              Netflix and <a href="https://engineeringblog.yelp.com/2014/10/using-moe-the-metric-optimization-engine-to-optimize-an-ab-testing-experiment-framework.html">Yelp</a> use Metrics Optimization software like <a href="http://github.com/Yelp/MOE">Metrics Optimization Engine (MOE)</a> which take advantage of Parallel Bayesian Optimization<d-cite key="yelpBO"></d-cite>.
            </p>
          </li>
        </ul>
      </li>
    </ol>

    <d-bibliography src="references.bib"></d-bibliography>
  </d-appendix>

</body>
</html>
